Distribution of Diffusion Constants and Stokes-Einstein Violation in supercooled 

liquids 
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It is widely believed that the breakdown of the Stokes-Einstein relation between the translational 
diffusivity and the shear viscosity in supercooled liquids is due to the development of dynamic 
heterogeneity i.e. the presence of both slow and fast moving particles in the system. In this 
study we directly calculate the distribution of the diffusion constant for a model system for different 
temperatures in the supercooled regime. We find that with decreasing temperature, the distribution 
evolves from Gaussian to bimodal indicating that on the time scale of the a relaxation time, mobile 
(liquid like) and less mobile (solid like) particles in the system can be unambiguously identified. We 
also show that less mobile particles obey the Stokes-Einstein relation even in the supercooled regime 
and it is the mobile particles which show strong violation of the Stokes-Einstein relation. Finally, 
we show that the degree of violation of the Stokes-Einstein relation can be tuned by introducing 
randomly pinned particles in the system. 



INTRODUCTION 

In normal liquids, the shear viscosity (77) is related to 
the translational diffusion constant (D) of a particle dif- 
fusing through it via the Stokes-Einstein (SE) relation 
as D = £jp where c is a constant which depends on 
the details of the particle and boundary conditions and 
T is the temperature. Although, originally derived for 
a macroscopic probe particle in the hydrodynamic limit, 
the SE relation also holds for the self diffusion of liq- 
uid particles at high temperatures However, it has 
been shown extensively (5l-[20j| that when a liquid is su- 
percooled the SE relation breaks down i. e. the measured 
self diffusivity becomes much larger than the value pre- 
dicted by the SE relation. Often the shear viscosity is 
substituted by the relaxation time (r) which is cheaper 
to compute and more commonly treated in theories. 

Phenomenological arguments e.g. by Stillinger and 
Hodgdon [4] and by Tarjus and Kivelson 

[HI 

have shown 

that by considering supercooled liquids to consist of mo- 
bile "fluid-like" and less mobile "solid-like" regions, the 
decoupling between the translational diffusion (D) and 
the relaxation time (r) is explained naturally because 
the average diffusion constant is predominantly deter- 
mined by the "fluid like" regions whereas the average 
relaxation time is dominated by the "solid-like" regions. 
The existence of clusters (with finite lifetime) of mobile 
and less mobile particles has been directly shown in many 
different studies [2T '2'A and is known as the dynamical 
heterogeneity (DH). However, the above-mentioned phe- 
nomenological theory of decoupling is based on the ex- 
istence of a distribution of diffusivity (relaxation times) 
whereas previous studies on DH have typically analyzed 
the distribution of particle displacements. Hence not 
much direct information about the existence and nature 
of the diffusivity (relaxation time) distribution is avail- 
able. 



There are however, indirect evidences, e.g. the univer- 
sal exponential tail in the Van Hove functions [24| seen 
for many supercooled liquids, which supports the exis- 
tence of a distribution of diffusivity (relaxation times). 
Consider an extreme case where the system has regions 
with two diffusivity - one for "solid like" (D s ) and the 
other for "liquid like" regions (D{). Hence a distribution 
of diffusivity can be written as p(D) = AS(D S ) + B5(Di) 
where A and B are fixed by the normalization condition 
and the amount of solid like and liquid like regions. Now 
if we calculate the van Hove correlation function as 



G s (x,t) = J dDp(D)g(x\D), 



(1) 



where g(x\D) = % / 4 ^ rDt ex P f — ISt) then one may show 
that the van Hove function will have a long tail and de- 
pending on the distribution of the p{D), the tail of the 
distribution can be either exponential or Gaussian [251 . 
In general the exponential tail has been reported [24| 
which, as mentioned in 25], might be due to the small 
range of the data. 

The main aim of this study is to calculate directly the 
distribution of the diffusivity from the simulation data. 
From the previous discussion, one may expect that in 
deeply supercooled liquids the distribution of diffusivity 
(p(p)) will be bimodal in general. Unlike previous works 
[26H28j where the presence of two types of particles with 
distinct mobilities are inferred indirectly from the dis- 
tribution of displacements, our method directly and un- 
ambiguously shows that the distribution of the diffusivity 
becomes bimodal below some temperature. Note that the 
bimodal nature of the distribution only says that there 
are two types of particles in the system, but it does not 
prove that these particles are clustered together to form 
"solid like" and "fluid-like" regions. Hence bimodal dis- 
tributions of diffusivity alone is not enough to justify the 
picture of supercooled liquids being a sparse mixture of 
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"fluid like" and "solid like" regions. Using random pin- 
ning geometry we further show that indeed one has "fluid 
like" and "solid like" regions in the system up to some 
time scale of the order of a relaxation time, r a . 

We end this section by noting that the phenomeno- 
logical explanation of decoupling does not tell anything 
about the nature and origin of the heterogeneity. Mi- 
croscopic theories to understand decoupling from the dy- 
namic heterogeneity e.g. mode coupling theory (MCT) 
I29L random first-order transition (RFOT) theory (30I — 
32l | . shear transformation zone (STZ) theory [33|, dynam- 
ical facilitation [34[ and the obstruction model |35| do not 
mutually agree on the origin and the nature of heteroge- 
neous dynamics. Besides, the decoupling of D and r 
may also be explained in terms of a growing length scale 
27 , [36 1 which does not directly require a distribution of 
diffusivity (relaxation times). 



DISTRIBUTION OF DIFFUSIVITY 



where p{D)dD = P(Dt)cL(Dt). The choice of Dt as 
our variable is due to the fact that D changes by several 
orders of magnitude in the studied temperature range 
whereas Dt changes relatively modestly with decreasing 
temperature and it will be easier to compare the distri- 
bution obtained for different temperatures. 




Below we explain briefly the method used to extract 
the distribution of diffusivity directly from the van Hove 
correlation function G s (x,t) using the iterative algo- 
rithm suggested in [37} and recently used in [25[ for the 
diffusion processes in biological systems. We start with 
the definition of van Hove correlation function 



G s (x,t) = (S [x - (Xi{r) - Xiffl))]) , 



(2) 



where the < . > implies averaging over the time origin. 
Now if we assume that particle displacements are caused 
by diffusion processes and there is a distribution of local 
diffusivity p(D), then we can express G s {x,t) in terms 
of p(D) as 



G s {x,t) 



D 



p(D).g(x\D).dD, 



(3) 



where g(x\D) = 1/V 4ttDt exp (— x 2 /ADt) and D is the 
upper limit of diffusion constant and will be equal to 
diffusivity for a free diffusion. Now given the G s (x,t) 
we calculate the distribution of diffusivity p(D) following 
37] as 



p n+1 (D)=p n (D) 



G s {x,t) 

— 00 7 ") 



g(x\D)dx, (4) 



where p n (D) is the estimate of p(D) in the n th iteration 
with p°(D) = (l/D a v g )exp(-D/D avg ) and 



/■A) 

G:(x,t)= / p n (D).g(x\D).dD. 
Jo 



(5) 



Similarly 



P n+1 (Dr) = P n (Dr) / 

J — I 



G s (x,t) 
G"(x,t) 



g(x\D)dx, (6) 



FIG. 1: The solid line shows the results of the iterative scheme 
used to calculate the distribution of diffusivity ( see text for 
details ) along with the original distribution shown as sym- 
bols. The agreement is really encouraging and inset shows 
the corresponding comparison for the van Hove correlation 
functions. 

We tested on a toy model whether the above men- 
tioned iterative scheme converges to the correct solution. 
We started with a distribution P(Dt) and calculated the 
van Hove correlation function G s (x,t) using EqUH and 
then used this to recalculate the probability distribution 
P(Dt) using Eq|6] In Fig|T] we have compared the ob- 
tained distribution P(Dt) with the input distribution. 
The agreement is really amazing with moderate number 
of iterations. In the inset we have compared the van 
Hove function obtained from the converged distribution 
of diffusivity with the exact one. Here also the agreement 
is near perfect. Note that the iterative scheme does not 
depend at all on the initial guess distribution. 

After establishing the rapid convergence of the itera- 
tive scheme to the correct solution we tried to calculate 
the distribution of diffusivity for a model glass forming 
liquid, the Kob- Andersen Model [38|. To compare the 
van Hove functions for different temperatures we cal- 
culated it for the time r m when the mean square dis- 
placement becomes half the inter particle diameter i.e. 
< Ar 2 (r m ) >= 0.50 as shown in FigJ5J The temperature 
dependence of r m is very different from the a-relaxation 
time T a calculated from the normalized two point corre- 
lation function as < Q(r a ) >— 1/e. In the inset of FigEJ 
we have shown r m and r Q as a function of temperature. 
This is another manifestation of the Stokes- Einstein (SE) 
violation in this model. 

To extract the distribution function of the diffusivity 
one needs to supply somewhat smoothly averaged data of 
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FIG. 2: Mean Squared displacement (MSD) for different tem- 
perature and the horizontal line indicates the time where 
mean squared displacement reaches 0.50 in the units of par- 
ticles diameter for different temperature. The corresponding 
time is denoted as r m here (see text for details). Inset: Com- 
parison of the a-relaxation time r a calculated from the nor- 
malized two point correlation function as < Q{T a ) >= f/e 
with r m . 



the diffusivity P(DT m ) as function of temperature and in 
Fig|4l we showed the van Hove correlation functions cal- 
culated from these distributions of diffusivity along with 
the simulation data. The agreement between the simula- 
tion data and the calculated ones is fantastic. An impor- 
tant point to mention here is that tails of these van Hove 
correlation functions can not be completely described by 
a single exponential function over the whole range at least 
for the low temperature data ( T < 0.50). Rather they 
can be better fitted by two exponential functions. 
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FIG. 3: Calculated distribution of the diffusivity for different 
temperatures. Notice the appearance of the bimodality in 
the distribution just below the onset temperature T = 1.00. 
Clear bimodal distribution in the supercooled regime confirms 
that there are two different types of particles in terms of their 
mobility up to time scale of typical relaxation time. The 
appearance of more peaks in the distribution at still further 
lower temperature is really interesting, indicating possibility 
of extremely slow to moderately slow to very fast particles in 
very deep supercooled state. 

the van Hove correlation function G s (x, r m ) for the iter- 
ative scheme to converge fast. In this study we fitted the 
extreme tail of the calculated G s (x,T m ) using exponen- 
tial function as this part is in general noisy and difficult 
to average. In Fig|3J we have shown the distributions of 



FIG. 4: The van Hove correlation functions (solid line) as 
obtained from the iterative scheme (see text for details) along 
with the simulation data (symbols). The curves are shifted 
upward for clarity and to point out that the line goes though 
the data points over the whole range of the data. 
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FIG. 5: The temperature dependence of the diffusivity asso- 
ciated with the solid like (D s ) and liquid like (Di) particles. 
These values are calculated from the peak positions of the 
distributions in Fig[3] Inset : The Stokes-Einstein violation 
parameter 6(T) = for the two types of particles. One 
sees clearly that solid like particles obey Stokes-Einstein re- 
lation to a reasonable accuracy over the whole temperature 
range and it is the liquid like particles which show strong 
Stokes-Einstein violation. 

Now looking at FigEl we can see that just below the 
onset temperature (T = 1.00), the distribution starts to 
become bimodal and two peaks clearly emerge at temper- 
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FIG. 6: Stokes-Einstein violation parameter has been shown 
as a function of temperature for different density of the frozen 
particles pi mp - One can clearly see that the deviation from 
the Stokes-Einstein relation becomes stronger with increasing 
density of the frozen particles. 



ature around T = 0.60. At further lower temperatures 
the distribution seems to show existence of shoulder or 
another peak but the peak at large diffusivity remains 
intact although with decreasing peak height. Thus we 
have clearly demonstrated that there are two types of 
particles in the supercooled liquid. One is more solid like 
(less mobile) than the other in the time scale of the or- 
der of a relaxation time, r a . It is worth noticing that 
the width of the distribution increases with decreasing 
temperature which indicates increase of DH leading to 
stronger SE breakdown [l7j. 

In the top panel of FigJSJ we have shown the temper- 
ature dependence of the diffusivity associated with the 
solid like (D s ) and the liquid like (£>/) particles. Here we 
took the peak position as an estimator for the diffusion 
constants for solid and liquid like particles. At the low- 
est three temperatures shoulders seem to appear in the 
distributions, but we calculated the peak position to be 
the position of the dominant peak. In the lower panel of 
Fig|5j we calculated the Stokes-Einstein violation param- 
eter 9(T) = for the two set of particles. One clearly 
sees that solid like particles obey the Stokes-Einstein re- 
lation over the whole temperature range, whereas the 
liquid like particles show strong SE violation leading to 
overall violation of the SE relation in the liquid. A simi- 
lar conclusion is reached in [26] . 

After establishing the fact that two types of particles 
can be identified in the system on a time scale of a re- 
laxation time, we now turn to the question of whether 
the "solid-like" particles group together to form clusters. 
To answer this question we performed another set of sim- 
ulation experiments where we randomly freeze (random 
pinning geometry) some of the particles and then ask 
what is the effect of this on the dynamics of the system. 
If the "solid like" particles are clustered together to form 
a region then if we randomly freeze some of the parti- 



cles in the system, there will always be instances where 
these frozen particles are part of these "solid like" regions. 
In that case, due to the frozen particles the relaxation 
of these regions will be hindered further and the relax- 
ation time of the whole system will increase dramatically 
with increasing density of these frozen particles. How- 
ever, these frozen particles will have very small effect on 
diffusivity which is mainly governed by the "fluid like" 
particles, so the diffusivity will not change that dramat- 
ically. In this scenario one expects to see an enhance- 
ment of the Stokes-Einstein breakdown if one increases 
the number density p imp of the randomly frozen particles. 
On the contrary if the clusters of "solid-like" particles are 
not present in the system then we expect that changing 
the density of randomly frozen particles will have little 
effect on the relaxation dynamics of these regions i.e. 
no dependence of the SE violation parameter with pi mp . 
In FigHJ we have shown the temperature dependence of 
the Stokes-Einstein violation parameter 9(T) = for 
different pinning densities pi mp . The deviation of 9(T) 
from being a constant as a function of temperature with 
increasing pi mp is very dramatic. 
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FIG. 7: Diffusivity (D) plotted as a function of relaxation 
time (r a ) for different density of the frozen particles pt mp 
in log-log to show the power law relationship between these 
quantities expected from the fractional Stokes-Einstein Re- 
lation D oc t c 7 1+ ", with lo > 0. Notice that the exponent 
uj increase with increasing pinning density pi mp indicating 
a stronger breakdown of the Stokes-Einstein Relation. The 
solid line has a slope equal to -f .0. 

In literature people often use a fractional Stokes- 
Einstein relation 0, [3 35 1 to fit the data where the 
original relation breaks down. We tried similar fits to 
the low temperature data using the form D oc t q T 1+ " ; , 
with u) > 0. In Fig 13 we have shown that our data can 
be well represented by the fractional Stokes-Einstein rela- 
tion in the low temperature range with u increasing with 
increasing pi mp . Now it is important to mention that 
in the past many people have tried to tune the Stokes- 
Einstein breakdown exponent to better understand the 
physics behind it. In those studies either different inter- 
action potentials [39( have been used or different spatial 
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dimensions was used 40-42] to see the change in the ex- 
ponent. The motivation for going to higher dimensions is 
that in the mean field limits one expects to have w = 0. 
Our study provides an new and easy way to tune this 
exponent. 

To conclude, we have shown by explicitly calculating 
the distributions of diffusivity for different temperatures 
for a model glass former in the supercooled regime that 
the state of the system can be well described by a sparse 
mixture of "fluid like" regions in the matrix of "solid 
like" regions on the time scale of a relaxation time. We 
also showed that the "solid like" regions to a great ex- 
tent follows the Stokes-Einstein Relation over the whole 
temperature range. All the deviation comes from the dy- 
namics of the "fluid like" particles leading to a overall 
breakdown of the relation. Finally we showed how with 
random pinning one can drastically enhance the decou- 
pling between the translational diffusion and the relax- 
ation time, thereby providing a new and easy way to tune 
the exponent associated with it to understand the phe- 
nomena even better. It will be interesting to extract the 
length scale associated with this "solid like" and "liquid 
like" region and compare that with the dynamic hetero- 
geneity length scale and the static length scales [HI , work 
in this direction is in progress and will be published else- 
where. 

We want to thank Srikanth Sastry and Chandan Das- 
gupta for many useful discussions. 
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